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We study the behaviour of scalar perturbations in the radiation-dominated era of Randall- 
Sundrum braneworld cosmology by numerically solving the coupled bulk and brane master wave 
equations. We find that density perturbations with wavelengths less than a critical value (set by 
the bulk curvature length) are amplified during horizon re-entry. This means that the radiation 
era matter power spectrum will be at least an order of magnitude larger than the predictions of 
general relativity (GR) on small scales. Conversely, we exphcitly confirm from simulations that the 
spectrum is identical to GR on large scales. Although this magnification is not relevant for the 
, cosmic microwave background or measurements of large scale structure, it will have some bearing 

on the formation of primordial black holes in Randall-Sundrum models. 
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I. INTRODUCTION 

Our view of cosmology has been revolutionized by the notion that the universe may be a lower-dimensional ob- 
. ject embedded in some higher-dimensional space. A particularly simple realization of this idea is furnished by the 
Oh' Randall-Sundrum (RS) braneworld model fT], which postulates that our observable universe is a thin 4-dimensional 
Q , hypersurface residing in 5-dimensional anti-de Sitter (AdS) space (see Refs. @, |3, 01 for reviews). Ordinary matter 
• degrees of freedom are assumed to be confined to the brane, while gravitational degrees of freedom are allowed to 
c/3 propagate in the full 5-dimensional bulk. The warping of AdS space allows us to recover ordinary general relativity 
^ (GR) at distances greater than the curvature radius of the bulk t. Current laboratory tests of Newton's law constrain 

[ £ to be less than around 0.1mm [5]. 
—{ . The cosmological implications of the Randall-Sundrum scenario have been extensively studied. It is well known that 
^ ' the Friedmann equation governing the expansion of the brane universe differs from general relativity by a correction 
'/"J , of order p/cr, where p is the density of brane matter and a > (TeV)^ is the brane tension. The magnitude of this 
correction defines the "high-energy" regime of braneworld cosmology as the era when p > a or equivalently H£ > 1 , 
where H is the Hubble parameter. At high energies, the RS Friedmann equation implies H oc p, which results in 
dynamics significantly different from the standard cx p expansion law. 
1^ I The major outstanding issue in RS cosmology is the behaviour of cosmological perturbations @, 0, The 
equations of motion governing fluctuations of the model are found to differ from GR in two principal ways at early 
times: First, they acquire 0{p/(j) high-energy corrections similar to those found in the Friedmann equation. By them- 
selves, such corrections are not difficult to deal with: they just modify the second-order ordinary differential equations 
(ODEs) governing perturbations in GR. But the second type of modification is more problematic: Perturbations on 
the brane are also coupled to fluctuations of the 5-dimensional bulk geometry, which are collectively known as the 
■ "Kaluza Klein" (KK) degrees of freedom of the model. The KK modes are governed by master partial differential 
' equations (PDEs) defined throughout the AdS bulk [1, 01 ■ In the original work of Randall & Sundrum [l| , the brane 
had a simple Minkowski geometry and the KK mode master wave equations were solvable via separation of variables. 
But the motion of the brane in the cosmological case breaks the time-translation symmetry of the bulk, which makes 
a simple separable solution unattainable in most cases. The exception is the case of a brane undergoing de Sitter 
inflation, where many analytic and semi-analytic results are now available pll. IT3. [isl. [fl . 

Many authors have considered various schemes to solve the perturbation problem without dealing with the KK 
master PDEs directly. The most straightforward approach is to simply set the KK degrees of freedom to zero 
while retaining ©(p/cr) corrections, which we refer to as the "4-dimensional effective theory". There are several 
approximation methods that attempt to move beyond this simplest effective theory in the case of tensor and scalar 
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perturbations [H, [13, [H, [H, H HH, [H [H, H, [H [H, [l^ . These generally consider the behaviour of fluctuations 
during certain limiting regimes (i.e, high-energy or low-energy) , but the only known way of tackling the problem on 
all scales simultaneously is by direct numerical solution of the equations of motion. 

The problem of tensor perturbations is somewhat simpler than the scalar case. The reason is that in linear 
theory, tensor modes are pure gravitational degrees of freedom that do not couple to matter fluctuations (provided 
that the matter anisotopic stress is neglected). In the braneworld context, this means that there are no brane- 
confined tensor degrees of freedom. The perturbation problem formally reduces to solving a wave equation in the 
bulk with boundar y co nditions imposed on a moving brane. Numeric solutions have been obtained by several groups 
m, IHj Hdl HH, [H,]!!, m, HH. interestingly, the effective theory predicts a blue spectrum (i.e. excess power) in the 
stochastic gravitational wave background at frequencies corresponding to modes that enter the Hubble horizon in 
the high-energy regime. However, once KK effects are incorporated via numeric simulations, one finds that the true 
spectrum is flat and virtually identical to GR. That is, the magnitudes of the 0{p/(j) and KK effects are nearly the 
same, but they act in opposite sense and end up canceling each other out. Physically, the reason that the KK modes 
cause a suppression of the amplitude at high frequencies is that they efficiently radiate gravitational wave energy from 
the brane into the bulk. 

The purpose of the paper is to numerically solve for the behaviour of scalar perturbations in the radiation-dominated 
regime of braneworld cosmology. (Numeric analysis of the scalar perturbations during inflation has already been done 
in Refs. [H, [3].) Unlike the tensor case, there are several scalar degrees of freedom residing on the brane, such as the 
density contrast 5p/p. The problem reduces to the solution of a bulk wave equation coupled to a master boundary 
field satisfying its own second order ODE on the moving brane. Recently, two of the numerical codes developed to 
deal with tensor perturbations [H, Isl] have been generalized to handle boundary degrees of freedom [3, |3a| ■ We 
use both codes in this paper, which gives us the ability to confirm the consistency of our numeric results via two 
independent algorithms. We are ultimately interested in finding the matter transfer function in the radiation era, and 
also determining the relative influence of KK and high-energy effects on the density perturbations. Heuristically, we 
may expect the KK modes to amplify high-energy/small-scale density perturbations, which is the opposite effect from 
the tensor case. The reason is that we know that the gravitational force of attraction in the RS model is stronger than 
in GR on scales less than I [l], [s^l • This implies that modes with a physical wavelength smaller than I during horizon 
crossing will be amplified due to the KK enhancement of the gravitational force. However, this physical reasoning 
needs to be tested with numeric simulations. 

The layout of the paper is as follows: The background geometry of the RS cosmology we consider is described in 
SjlTl The formulae governing the gauge invariant scalar perturbations are given in i jllll Also in that section, we derive 
new analytic approximations for various perturbation variables at very high energies. In ijlVl we present the numeric 
algorithms that we use to solve the coupled system of bulk and brane wave equations. Our results for the matter 
transfer function are also contained in that section. Finally, we summarize and discuss the implications of our work 
in gV] 

II. BACKGROUND GEOMETRY AND BRANE DYNAMICS 

As discussed in the Introduction, the Randall-Sundrum model we will be considering consists of a thin brane, which 
is realized as singular 4-surface, residing in 5-dimensional anti-de Sitter space. In Poincare coordinates, the bulk 
metric is given by the line element 

dsj = — (-dr^ + Sijdx'dx^ + dz^). (1) 

Here, i is the curvature scale of the bulk. The boundary of this spacetime given by a brane is described parametrically 
by 

r = rb(77), z = Zb(77). (2) 

The parameter ry is defined to be the conformal time along the brane, which leads to the following induced line element 

d-q^ ^ drl ~ dzl, dsl ^ a^{-drf + dx^) = -dt'^ + a^dx^ . (3) 

Here, we have identified the scale factor as 0(77) = 1/ Zh{rj) and t as the cosmic time. Now, we define m° — dx^/dt to 
be the 5-velocity field of comoving observers on the brane, while n" is the brane normal. These vector fields can then 
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be used to construct directional derivatives parallel and orthogonal to the brane: 

d „.d\ d Id 



du = u-da = - f + HH^^ - Hi-^] =^^-^, (4a) 
a \ OT oz J at a orj 

dn = n'^da = - {-Hl^ + VTTW|-) , (4b) 

a \ OT oz J 

respectively. Here, H is the usual Hubble parameter 

fdur\ Ida r £ 

where we have introduced the alternative bulk coordinate r. We assume that the matter on the brane has a fluid-like 
stress-energy tensor with density p and pressure p: 

Tab = {P + p)UaUb + phab, Qab ^ Kb + naUb- (6) 

The junction conditions are given in terms of the extrinsic curvature Kab of the brane, 

Kab^ha'V.nb, 2Kab = -4[Tab-^KbiT-a)], (7) 
where Kg — is the gravity-matter coupling in 5 dimensions and a is the brane tension. This yields 



= ■ i^) 
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We assume the standard Randall-Sundrum fine-tuning 

kI = kII^SttGI^ (9) 

to obtain the Fricdmann equation 

1 p / p\ 8ttG 



Note the high-energy 0{p/a) correction to the expansion rate. Finally, note that we have the usual conservation of 
stress-energy on the brane; i.e., 

/i^'VaTbc = 0, ^ = -3(1 + u.)pH, w = ^. (11) 
at p 



III. SCALAR PERTURBATIONS 
A. Gauge invariant built perturbations 

We now turn our attention to the perturbation of the cosmology introduced in It has been shown in Refs. [1, 0| 
that scalar-type perturbations of the bulk geometry |T]) are governed by a single gauge invariant master variable Q.. 
We summarize those results in this subsection. 

First, we introduce a harmonic basis with mode functions^ 

y = e^'^■^ = Y,, ^ ^d^d.Y +^5,jY. (12) 



^ To describe the perturbations in the RS model we mostly follow the notation and conventions of Kodama et al. 0, which are based on 
the 4-dimensional formalism of Kodama and Sasaki .38| . 
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Using these mode functions, we write the perturbed bulk metric as 

„2 



(13) 



where 



Qafjdx^dx^ = — (-dr^ + dz"^ 



(14) 



Implicit integration over the wavevector k is understood in equations like (|T3|) . A priori, the particular form of the 
perturbation variables faf3, fa, and Ht depend on the choice of gauge. However, one can construct gauge invariant 
combinations as follows 



O T 



where 



(15) 



(16) 



and Da is the covariant derivative associated with the 2-metric qa(j- The gauge invariant quantities are then given in 
terms of the master variable — ^l{t, z) 



F 



1 

6r 



The bulk master variable satisfies the following wave equation 







d^n 3 dn 






dz^ z dz 


(z2 



(17) 



(18) 



B. Gauge invariant brane perturbations 



In addition to perturbations of the bulk metric, we must also consider perturbations of the brane metric. These 
are parameterized as 

dsl = (1 - '2aY)dt^ - 2a(3Y,dt dx' + [(1 + 2/iL)(5y + 2hTY,j] dx'dx'. (19) 
As in 4 dimensions, the gauge-dependant variables (a, /?, h^, hx) can be used to construct gauge-invariant quantities 

1 Ha 
3 

where 



(20) 



a dhr 



(21) 



As demonstrated in Ref. [s^ , the effective Einstein equation on the brane is in general given by 



n„ 



1, 



-T 



1 



rTT„, 



1 



Tcd.T'^'^hab 



-T^h, 



(22) 



1 

„ . -^2""" ' 8"""~ 24^ 
Here, £ab represents the projection of the electric part of the bulk Weyl tensor onto the brane. Notice that £ab satisfies 

h'^^'Eab = 0, Va(«;iH"'' - = 0, (23) 



and vanishes in the background geometry discussed in gH] When perturbing (j22|) . one can treat Sab as an additional 
fluid source with a radiation-like equation of state. Hence, we parameterize the perturbations of Tab and Sab as 



5T,^ =6pY6,^ +k^6TTY,^, 



6£o° = 46p£Y, 
6£° = KlkYiSqs, 

5£^ = -t4{\5p£Y5^^ + kH^£Y^^). 



(24a) 
(24b) 
(24c) 
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Henceforth, we will assume that the matter anisotropic stress vanishes Sn = 0. The Weyl fluid perturbations 
{Sp£,6q£,STT£) are the "Kaluza-Klein" (KK) degrees of freedom of the model alluded in the Introduction, since 
they represent the efl'ects of bulk geometry fluctuations on the brane. The matter perturbation variables above are 
not gauge invariant, but the following quantities are: 



pAg = Sp£ - 3HSq£, 



a{p + p)V 
a{p + p)V£ 



-kSq - {p + p)aag, 
-kSq£ - {p + p)aag. 



r^Sp- clSp, (25a) 
(25b) 



Here, we have defined the sound speed as — 5p/5p. Since the Weyl fluid has the equation of state 5p£ — \5p£, 
there is no KK entropy perturbation Fg = 0. Furthermore, the KK anisotropic stress is automatically gauge invariant. 
Armed with these definitions, it is possible to derive the following gauge- invariant equations from the perturbation of 
Eq. dill): 



$ = 



AirGpa^ 



1 



jj^ d^^ AT^G{p+p)a 
dt k 
$ + * = -STrGa^fef. 



A. 



V£ 



(26a) 

(26b) 
(26c) 



The first of these gives the high-energy 0{p/a) and KK corrections to the Poisson equation on the branc. The second 
equation gives the modifications of the standard equation governing the evolution of the peculiar velocity field. The 
third equation demonstrates how the KK-modes can give rise to an anisotropic stress on the brane. 

Additional equations can be obtained from perturbing the equation representing the conservation of matter stress 
energy on the brane: 5{h°'^V aTbc) = 0. These are 



{aV) 



a dt 

4-^(«VA) = --p(l 
a"^ dt a 



kV + clpA 
a ' a p(l + w) 



w) 



3a2 dH 
'dt 



U;) ( ^ _ i/^ 

^ * dt 



(27a) 
(27b) 



Finally, we can also perturb the righthand expression in Eq. (j23p to get two more equations, but these are not sufficient 
to close the effective Einstein equations (1^^ on the brane. Hence we can only go so far with this effective Weyl fluid 
description; a complete treatment must describe how the brane degrees of freedom presented here are coupled to the 
bulk degree of freedom Q. 



C. Perturbation of the junction conditions 



The means to connect the bulk perturbations of mil Al and the brane perturbations of ijHIBI is the perturbation of 
the junction conditions ([7]). This yields several results, the most important of which is that f2 satisfies a boundary 
condition on the brane 



1 



6pa^ 



A 



0. 



The perturbation of ([7]) also gives the KK gauge invariants in terms of ft: 



1 \ ^ 

-4PA. = 



Kla{p + p)V£ = 



fc3 



3£a2 
1 

2£^ 



3a2 



[Hn - du^] 



b ' 

Ik^n - Hdun - ^^ d^n + din 

3 icr 



(28) 

(29a) 

(29b) 

(29c) 
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These reproduce the results first derived in Ref. [40'|. Using these equations, we can find expUcit formulae for the 
gauge invariants defined above in terms of A and Qhiv) — ^(■''b('7)j -ZbC??)): 



3pa'^{3wp + 4:p + a) 



2£a^ drj ' 
(3w + 4)/92 (5 + 3w)p fc2 



3wHa ^ 
V = T- -A - 



2Pc 



1 dA kH 



2Paa 
k dQh 



3£a3 



3H df7b 



2£a^ drj 2£a^ dj]^ 



fc(l + w) k{l + w) dr] 2£a^ 2la^ drj ' 
Other quantities of interest are the curvature perturbation on uniform density slices, 



HaV 



A 



1 3pa?{wa — a — p) 
3 " " 



k 3(1 + w) 
and the curvature perturbation on comoving slices, 

HaV 'ipa^{wa ~ a — p) A 



A 



Ha dA 



l + w k'^{l + w) dr] 6£a^ 



Ha dA 



k^£^cr^ 



l + w fc2(l + w) drj 6£a^ 



(30a) 
(30b) 
(30c) 

(31) 

(32) 



Hence, if we have knowledge of the density contrast A and the value of on the brane, we can obtain the behaviour 
of all the other gauge invariants via algebra and differentiation. 

The above formulae can be easily manipulated to find a wave equation for A: 



d^A , o , dA 
-— + (1 + 3cl - 6w)Ha— 
drj'^ drj 



2,2 3pa^ ^ ip^a^ ^ 
cik^ + + „„„ B 



a£^ 



2p 



A = - 



p 



3£a^ 



^ = 6c2-l-8w + 3w^ S = 3c2-9w-4. 



(33a) 
(33b) 



The above ODE, the bulk wave equation (fT8|) and the boundary condition (|28|) comprise a closed set of equations for 
A and fi. 

Before moving on, we have three remarks: First, note that in the low energy universe, we can neglect 0{p^ /a^) 
terms. Then, the A wave equation reduces to 



d^A 



dA 



(1 + 3c^ - 6w)Ha— + \clP + AnGpa'^iGcl - 1 ~ 8w + 3w^)] A « + 

drj p 



k^T k'^{l + w)nb 



3£a^ 



(34) 



Setting r^b = yields the standard 4-dimensional dynamical equation for A; hence, we recover GR at low energies. 
Second, we note that from the conservation of stress-energy on the brane (|27p . we must have 



dC 

drj 



kV HaT 
~ ^ p{l + w)' 



(35) 



We can easily verify that Eqns. ()30c|) . (I3ip and (l33|) imply that ([35|) is satisfied identically, which is a good consistency 
check of our formulae. Finally, the entropy perturbation can be simplified when the brane matter is either a perfect 
fiuid or a scalar field: 



r = (T - c2)pA, T 



cj., perfect fluid. 



1, scalar field. 

For the rest of this paper, we will be considering the perfect fluid case. 



(36) 



D. Asymptotic behaviour in the high-energy radiation-dominated regime 

To model the perturbations in the very early universe, we assume radiation domination and p 3> cr; or, equivalently, 
H£ ^ 1. Under these circumstances, the following approximations hold: 



w = 1/3, a~ao{krjy^^, pftiaH£, 9„ 



-a, = -a 



(37) 
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The last relationship can be used to re-cast the boundary condition on the bulk master variable Q into a 1st order 
ordinary differential equation for fib- Hence, the equation of motion for A and boundary condition become 



By differentiating the first equation we can derive a decoupled third order equation for A, which can be solved via 
power series methods. The solutions for A and fl are a superposition of three linearly independent modes: 

3 3 
i=l 1=1 

where the Ai are constants. At leading order in (fcr;), these mode functions are given by 

dominant growing mode: A^^^ « |(^'?)^, ^ ~ a^k^'^ i{kr]f , (40a) 

subdominant growing mode: A^^^ « (fc??)*^^ ^^^^ ~ -|a^^fc"2^(fc77)l/^ (40b) 

decaying mode: A^^^ « ^Hag \fc?7)"\ f^jf^ w -20agfc"^€^ (40c) 

As the labels on the left suggest, the growth of density perturbations for the dominant growing mode is faster than 
for the subdominant growing mode on super horizon scales. The density contrast of the third mode decays in the 
high-energy regime. 

We can calculate the behaviour of C and ^> for the dominant growing mode by substituting the full series solutions 
for A(i) and f^^^) into ((30b| and (EJ), which yields: 

C(i)«i, ^^'^^~4C^'K (41) 

We see that the dominant growing mode curvature perturbation and Newtonian potential are conserved on super- 
horizon scales. Also note that ([4T|) implies 

A(i) = ±fAy^(i). (42) 

Eqs. (j4ip and (|42p are similar to the standard super horizon growing mode results for a radiation dominated universe 
in general relativity, 

,T,(i)^. ^(1)^. 3(l-fu;) 2(1) (1) ^ {3w + l)\l + w) ( k Vai) _ ^ ( k \\(,) 

but the numerical factors are different.^ 

Finally, we note that the method described in the subsection fails to yield an approximate solution for the other 
metric potential <i>(^^ The reason is that when the full series expansions for A(i) and are substituted into 
pOap . the leading order contribution vanishes, leaving a result that is the same order as the error in the original 
approximations p7p . We must therefore rely on numeric simulations to determine the asymptotic behaviour of $. 

IV. NUMERIC ANALYSIS 

A. Dimensionless parameters and integration algorithms 

Our goal in the section is the solution of the system comprised of the bulk wave equation (|18[) and the brane wave 
equation (|33p subject to the boundary condition (pS)) . For the rest of the paper, we will restrict ourselves to the case of 



^ It is interesting to note that it is impossible to find an effective equation of state parameter uictt to make the general GR formulae 
given in II43I I compatible with II41I I and 11421 . Some authors |41| have previously tried to describe the high-energy radiation epoch of RS 
cosmology with an effective fluid with uictt = 5/3, but we see that this would predict A'^) 2i /b(k / H a)'^ C,^^^ and l'(l) -4C(i)/5, 
which is in clear conflict with the correct results l|41^ and II42II . This is due to a large Weyl anisotropic stress that modifies the GR 
relationship between ^ and <J>. 
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a radiation-dominated brane with w = 1/3. To perform the analysis, it is necessary to make the various quantities in 
the equations dimensionless. To do so, we define the "*" epoch as the moment in time when a mode with wavenumber 
k enters the Hubble horizon 

k ^ H^,a^:, i/* = H{r]^), a, — a{T]^), — Zb(?7*)- (44) 

Then, we introduce dimensionless/normalized quantities as 

n = p=^, H = He, k = H, = — = kz, = V/5*(/5* +2), a - — . (45) 

Another important era is the critical epoch when He = Hc( = 1 and the radiation density has its critical value 
Pc — — 1. We define the critical epoch as the transition between high and low energy regimes. The ratio of the 
wavenumber of any given mode to the wavenumber k^ ~ H^ac of the critical mode that enters the horizon at the 
critical epoch is 

A = Vp*(p* + 2) ' (46) 

where Oc is the value of the scale factor at the critical epoch. Generally speaking, we call modes with k > kc 
"supercritical" and modes with k < k^ "subcritical" . The scale defined by the critical mode in today's universe (with 
scale factor ag) is given by 



= 1.4 X 10^" m, (47) 

where gc is the number of relativistic degrees of freedom in the matter sector at the critical epoch. For £ — 0.1 mm, this 
corresponds to a scale of ^ 10 astronomical units (AU), which is incredibly tiny by cosmological standards. Finally, 
if we normalize all coordinates by z*, 

f = r/z*, z = z/z^, T) = -q/z^, (48) 

we find that the entire system of master equations is characterized by the dimensionless matter density at horizon 
crossing p*. 

Once the system has been reduced into the dimensionless form, we have two independent numerical codes that can 
be used to solve for A and f2b. The first is the pseudo-spectral (PS ) method used in Ref. [13, [s^l and the second 
is the characteristic integration (CI) algorithm developed in Ref. |3g]^; detailed descriptions of each method can be 
found in the respective papers. There is one technical difference between the two algorithms that is worth highlighting 
here: Namely, each code solves for VL in different regions of the 5-dimensional spacetime, as shown in Fig. [TJ One 
implication of this is that the initial data for each code needs to be specified at different places. For the PS code, one 
needs to choose Cl and dCl/dr on an initial spacelike hypersurface dM.~^ of constant f, while for the CI algorithm one 
needs Cl only on an initial null hypersurface SA^^i intersecting the brane. 

What initial conditions should we actually use? The calculations of §111 Dl suggest that at sufficiently early times, 
the dynamics of the system are well described by three distinct modes (j40p . If our initial data surface is positioned 
within this early era, it follows that generic choices of initial conditions will excite some superposition of these three 
modes. However, after a short period of time the dominant growing mode will overtake the contributions from the 
subdominant growing and decaying modes. Hence, we expect that the late time dynamics of our model will be 
insensitive to the choice of initial conditions, provided that the initial data hypersurface is oriented at an early enough 
time. Since we do not expect the initial conditions to matter very much, we make the simple choices 

PS initial conditions: A(0) = TVaf , 9^A(0) = 6JVH^aJ, n{dM^s) = 0' 5^r2(9Al;J = ejVp^a^pH,, (49a) 
CI initial conditions: A(0) =7Vaf, dr,A{0) = 6AAiJ,aJ, Cl{dMc^) = 0, (49b) 



^ To apply the CI method as described in [36l |. the bulk wave equation needs to be mapped into a canonical form via the change of variable 

n = {z,/zf/^<f,. 
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r 



PS computational 
domain 



simulation start 



>2 



FIG. 1; The computational domains employed by the pseudo-spectral (PS) and characteristic integration (CI) methods. The 
future and past boundaries dMp^ of the PS domain A4ps are lines of constant r, while for the CI method they are null rays. 
In the PS method, one needs a constant z regulator brane to render the computational domain finite. In principle, the position 
of the regulator is free, but it should be placed to the right of the point P in order to be outside the causal past of the physical 
brane. Initial conditions for the two methods are placed on dA4p^ and respectively. 



where d^j — d/df), — d/dr, and an "z" subscript denotes the initial value of the scale factor and Hubble parameter. 
Here, A/" is a normalization constant that we will often select to make C = 1 at the initial time. These initial conditions 
are motivated by the fact that we expect A oc 3> fib for the dominant growing mode at very early times, which 
can certainly be satisfied by setting f2 = on the initial data surface. For the PS method, the initial time derivative 
of Cl is selected to satisfy the boundary condition ([28]) at the initial time. Note that the initial conditions for the two 
methods are actually incompatible due to the different locations of the initial surface. But as we have argued above 
and will see below, this difference should make no difference to the final results. (We will test this assumption in 



We now turn to the results of our simulations. In Fig. [2] we show the output of the PS and CI codes for a 
typical simulation of a supercritical mode with p* =50. As expected we have excellent agreement between the two 
codes, despite the fact that they use different initial conditions In fact, we have confirmed that the codes give 

essentially identical results for a wide range of parameters, which gives us confidence in our numerical methods. It is 
also reassuring that the simulation results closely match the analytic predictions of ijlH Di on high-energy/super horizon 
scales. Note that for all simulations, we recover that A and Q are phase-locked plane waves. 



at sufficiently late times k-q ^ 1, which is actually the same behaviour as seen in GR, where we have the following 
exact solutions for the growing mode density contrast and curvature perturbation during radiation domination; 



flVD] below.) 



B. Typical waveforms 



A(77) (X cos ^p, A(?7) w 4((?7) 



(50) 




k-q < 1, 
fc?7 > 1, 



kr] < 1, 
k-q > 1, 



(51b) 



(51a) 
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FIG. 2: Comparison between typical results of the PS and CI codes for various brane quantities {left); and the typical behaviour 
of the bulk master variable {right) as calculated by the CI method. Very good agreement between the two different numerical 
schemes is seen in the left panel, despite the fact that they use different initial conditions. We also see excellent consistency 
between the simulations and the approximations developed in ijlll DI for the behaviour of the system on high-energy/superhorizon 
scales [c.f. Eqs. (|4Up - (|42|l ]. Also note that on subhorizon scales, A and C, undergo simple harmonic oscillations, which is consistent 
with the behaviour in GR (|5ip . The bulk profile demonstrates our choice of initial conditions: We see that the bulk master 
variable Cl is essentially zero during the early stages of the simulation, and only becomes "large" when the mode crosses the 
horizon. 



where ^ is a constant. 

In Fig. [21 we show the simulated behaviour of several different gauge invariants on the brane for given subcritical 
(fc — 0.25fcc) and supercritical (fc = lOfcc) modes. These plots illustrate several points that are generic to all of our 
simulations: The curvature perturbations TZc and C are conserved on superhorizon scales for all cases, just as in general 
relativity. This useful fact allows us to define the primordial amplitude of a given perturbation as the value of C or TZc 
before horizon crossing. For the metric perturbations $ and ^, wc find that for supercritical and superhorizon scales 

l^f] > l^l, -4C = constant, (for a < ac and fc < i?a). (52) 

This implies that the Weyl anisotropic stress k\5t:£ « —^/a^ is relatively large in the high-energy superhorizon 
regime. 

Finally, Fig. [5] illustrates how the ordinary superhorizon behaviour of perturbations in GR is recovered for modes 
entering the Hubble horizon in the low energy era. We see how A, $ and 4" smoothly interpolate between the non- 
standard high-energy behaviour to the usual expectations given by (|43[) . Also shown in this plot is the behaviour of 
the KK anisotropic stress, which steadily decays throughout the simulation. This is characteristic of all the cases we 
have studied. 



C. Enhancement factors and the transfer function 



If we examine the curvature perturbations plots in Fig. [3] in detail, we see that the amplitude of C, increases during 
horizon crossing. Furthermore, the degree of enhancement seems to increase with increasing k. This is quite different 
from the behaviour of the perturbations of a radiation dominated universe in ordinary GR (|5ip , where the asymptotic 
amplitude of C, is the same before and after horizon crossing. Hence, in the braneworld case we see an enhancement 
in the amplitude of perturbations that is not present in conventional theory. 

What is responsible for this enhancement? As discussed in the Introduction, there are actually two separate effects 
to consider: First, there is the modification of the universe's expansion at high energies and the 0{p/a) corrections to 
the perturbative equations of motion ([50)1 - ([55)1 . Second, there is the effect of the bulk degrees of freedom encapsulated 
by the bulk master variable f2 (or, equivalently, the KK fluid variables Ag, Vg and 5tts)- To separate out the two 
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FIG. 3: Typical behaviour of various brane gauge invariants for modes entering the horizon in subcritical {upper panels) and 
supercritical {lower panels) epochs. For the plots on the left, the grey straight lines give the analytic expectations (|40p for the 
behaviour of A and fib on large scales and high energies. The central plots show the comoving and uniform density curvature 
perturbations, TZc and respectively. As in conventional cosmology, these are conserved and equal on superhorizon scales. 
Also note that the amplitude of ( after horizon crossing is larger for the mode with larger k. The plots on the right show the 
two gauge invariant metric perturbations on the brane. Before the critical epoch we have that l^*] 2> |$|, while afterwards we 
recover the GR result $ ~ — vp. Finally, in the high-energy regime we see 4* « —4^, in line with the approximation (|4ip. 



effects, it is useful to introduce the 4-dimensional effective theory discussed above where all 0{p/cr) corrections to 
GR are retained, but the bulk effects are removed by artificially setting = 0. In the case of radiation domination, 
we obtain equations for the effective theory density contrast A^ff and curvature perturbation C^f,- from pip and (j33p 
with = 0: 

" U ^ 2^fc^ + J + ^^^^^ 

These in conjunction with the Friedmann equation pU]) and the conservation of stress energy pT|) give a closed set of 
ODEs on the brane that describe all of the 0{p/a) corrections to GR. 

In Fig. O we plot the predictions of GR, the effective theory, and the full 5-dimensional simulations for the behaviour 
of C and A for a supercritical mode. Since in any given model we expect the primordial value of the curvature 
perturbation to be fixed by inflation, it makes physical sense to normalize the waveforms from each theory such that 

CsD « Ccft ~ Cgr ~ 1, a<a,. (54) 

When this is enforced we see that the effective theory predicts a larger final amplitude for the density perturbation 
than GR. Furthermore, the final amplitude in the 5-dimensional simulation is larger than both of the other theories. 
From this we can infer that both 0{p/a) and KK effects induce enhancement in the amplitude of perturbations. 
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FIG. 4: The simulated behaviour of an extremely subcritical k <^ kc mode on superhorizon scales. On the left we show how 
the A gauge invariant switches from the high-energy behaviour predicted in Eq. (|42|l to the familiar GR result (143 |l as the 
universe expands through the critical epoch. We also show how the KK anisotropic stress n^Sive steadily decays throughout 
the simulation, which is typical of all the cases we have investigated. On the right, we show the metric perturbations $ and 'I' 
as well as the curvature perturbation (. Again, note how the GR result $ « —'i/ ~ —2(^/3 is recovered at low energy. 



This is in contrast to the situation for tensor perturbations in the Randall- Sundrum model, where modification 
of the expansion serves to increase the amplitude of gravitational waves while the bulk effects tend to decrease it 
[2i,[2l[3i|,[33,[3l,|3i,[3l. 

The amount of enhancement should be a function of the frequency of the mode, since we expect extremely subcritical 
modes k to behave as in GR. To test this, we can define a set of "enhancement factors", which are functions 

of k that describe the relative amplitudes of A after horizon crossing in the various theories. As in Fig. O let the 
final amplitudes of the density perturbation with wavenumber k be Csr>ik), C^ttik) and Ccn^fc) for the 5-dimensional, 
effective and GR theories, respectively, given that the normalization (jM]) holds. Then, we define enhancement factors 
as 

e^^'^liy- 

It follows that Qetf(fc) represents the 0{p/<j) enhancement to the density perturbation, Qeik) gives the magnification 
due to KK modes, while Q5D{k) gives the total 5-dimensional amplification over the GR case. These enhancement 
factors are shown in the left panel of Fig. [Sj We can see that they all increase as the scale is decreased, and that they 
all approach unity for k 0. Since Q = 1 implies no enhancement of the density perturbations over the standard 
result, this means we recover general relativity on large scales. For all wavenumbers we see Q^tt > Qs > 1: which 
implies that the amplitude magnification due to the 0{p/(j) corrections is always larger than that due to the KK 
modes. Interestingly, the Q-factors appear to approach asymptotically constant values for large k: 

Q.„(fc)»3.0, Q£(fc)w2.4, Q5D(fc)«7.1, fc » fee- (56) 

It is difficult to know if these are the true asymptotic limits for fc — s- oo due to the limitations of computing power. 

In cosmological perturbation theory, transfer functions are very important quantities. They allow one to transform 
the primordial spectrum of some quantity set during inflation into the spectrum of another quantity at a later time. In 
this sense, they are essentially the Fourier transform of the retarded Green's function for cosmological perturbations. 
There are many different transfer functions one can define, but for our case it is useful to consider a function T{k) that 
will tell us how the initial spectrum of curvature perturbations V'^^ maps onto the spectrum of density perturbations 
Va at some low energy epoch within the radiation era characterized by the conformal time rj > r/c- It is customary 
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FIG. 5: A comparison of the behaviour of the curvature perturbation (left) and the density perturbation A (right) in the 
fuU 5-dimensionaI theory including KK contributions (5D), the effective 4-dimension theory including 0{p/a) corrections (eff), 
and ordinary general relativity (GR). The waveforms for each theory are normalized such that = 1 on superhorizon scales. 
We can clearly see that given the same primordial values of the curvature perturbation, the final amplitude of the density 
perturbation in the 5D theory Car, is larger than in the effective theory Cotf, which in turn is larger than the GR value Cgr = 4. 




relative wavenumber log(/c/fcc) relative wavenumber log(/c/fcc) 



FIG. 6: Density perturbation enhancement factors (left) and transfer functions (right) from simulations, effective theory, and 
general relativity. All of the Q factors monotonically increase with k/kc, and we see that the A amplitude enhancement due 
to 0{p/a) effects Qctt is generally larger than the enhancement due to KK effects Qg. For asymptotically small scales k ^ kc, 
the enhancement seems to level off. The transfer functions in the right panel are evaluated at a given subcritical epoch in 
the radiation dominated era. The T functions show how, for a fixed primordial spectrum of curvature perturbations P^"', the 
effective theory predicts excess power in the A spectrum Va oc T^p'^"' on supercritical/subhorizon scales compared to the GR 
result. The excess small-scale power is even greater when KK modes are taken into account, as shown by Tar>{k; rf). 
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to normalize transfer functions such that T{k; rf) ^ > 1, which leads us to the following definition 



H{r])a{r]) 



Afc(r/) 



(57) 



Here, Qj^^ is the primordial value of the curvature perturbation and Akirj) is the maximum amplitude of the density 
perturbation in the epoch of interest. As demonstrated in Fig. [4l we know that we recover the GR result in the 
extreme small scale limit 



Afc(7?) 
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H{ri)a{r]) 



(58) 



which gives the transfer function the correct normalization. In GR, the transfer function is accurately given by 



1, k< 3Ha, 

[iHa/kf, k > 3Ha. 



(59) 



In the righthand panel of Fig. [51 we show the transfer functions derived from GR, the effective theory and the 5- 
dimensional simulations. As expected, the T(k; rj) for each formulation match one another on subcritical scales k < k^. 
However, on supercritical scales we have T^-^ > T^a > Tgr . 

Note that if we are interested in the transfer function at some arbitrary epoch in the low-energy radiation regime 
Ha 3> fee, it is approximately given in terms of the enhancement factor as follows: 



Ts,o{k\T]) 



1, 



k < 5Ha, 



pHa/kfQ,^{k), k>3Ha, 
Now, the spectrum of density fluctuations at any point in the radiation era is given by 



16 , / fc 



VA{k;v)^^T'ik;r^)l— r^\k). 



(60) 



(61) 



Using the enhancement factor results in Fig. [5] with Eqns. ([55)1 and ([50]) . we see that the RS matter power spectrum 
(evaluated in the low-energy regime) is ^ 50 times bigger than the GR prediction on supercritical scales k ^ lO'^/cc- 



D. Alternate initial data 



Before concluding, we wish to briefly revisit the issue of initial conditions. We have argued above that the existence 
of a dominant growing mode at high-energy and on superhorizon scales means that the details of initial data for our 
simulations are unimportant. To some degree, the fact that the PS and CI algorithms produce essentially identical 
results is a good confirmation of this, since each method uses a different prescription for initial data. 

But what if we were to use different classes of initial data? For example, we could modify the CI initial data such 
that the bulk field has constant amplitude along the initial data surface: 

(l{dMcj) = ^i = constant. (62) 

The waveforms generated by such a choice are shown in Fig. [T] Though the overall amplitude and early time behaviour 
of the signal seems to be sensitive to the initial value of VL, the ratio of the final amplitude of A to the initial value of 
C is the same for each choice of initial data. We have confirmed that this is also true for other choices of Cl{dM~^)^ 
including the case when the initial data is a simple sinusoid.'' Ultimately, it is the ratio /^[kr] ^ l)/({kri <C 1) that 
is relevant to the transfer function used to transcribe the predictions of inflation into observable quantities, hence we 
can be confident that our principal results hold for reasonable modifications of initial data. 



Sinusoidal initial data was used in Ref. |35| to test the insensitivity of the final spectrum of the stochastic background of gravitational 
waves to initial conditions in RS models. 
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FIG. 7: Simulated waveforms for Q = constant initial data in the CI method. Here, (f2/A)i is the ratio of the bulk field to the 
density perturbation at the initial time; this ratio is zero for the other simulations in this paper. The early time behaviour of A 
and Clh shows some sensitivity to {Cl/ A)i, but otherwise the waveforms are identical to one another up to an overall amplitude 
scaling. Indeed, if we normalize the simulation results by the initial value of we find that all of the A waveforms are coincident 
at late times, which means that the enhancement factor Qsd and transfer function T^u are insensitive to {Q,/A)i. Note that at 
early times we have oc a, which matches the behaviour of the subdominant growing mode from Eq. (|40|) . This confirms our 
expectations: varying the initial data causes the other subdominant modes to be excited to various degrees, but the dominant 
growing mode always "wins" at late times. 



V. DISCUSSION AND IMPLICATIONS 



In this paper, we have written down the equations of motion for generic gauge-invariant scalar cosmological pertur- 
bations in the RS braneworld model in a form suitable for numeric analysis We have developed new analytic 
approximations for the behaviour of fluctuations on high-energy superhorizon scales ( ijIIID[) for a radiation dominated 
brane. We have applied two different numerical algorithms to solve the equations in the early radiation-dominated 
universe f ijlVj) . Our numerical results show that the amplitude of modes which enter the Hubble horizon during the 
high-energy regime gets enhanced over the the standard GR result. Conversely, modes which enter at low energy do 
not show any late-time deviations from standard theory, as seen in Fig. 3) We only recover standard results for the 
metric perturbations, such as ^E" « — at low energies. 

Our simulations confirm the common sense expectation that all tangible effects from the fifth dimension are on 
scales smaller than a critical value k~^, whose physical size today is given by Eq. (|47p . To parameterize the degree 
of amplification as a function of scale, we introduced so-called "enhancement factors" and transfer functions in 
mV CI These show that the degree of enhancement over the GR case seems to reach a constant value at high k; the 
amplification of the subhorizon matter power spectrum is ^ 50 for k lO^fcc. We presented analytic arguments and 
numerical evidence that our results are robust against modifications of the initial data for simulations f ijlV Pp . 

As discussed in and i jlVCl the enhancement of perturbations upon horizon crossing can be attributed to two 
effects: namely, the 0{p/a) corrections to the background dynamics and the influence of the KK modes. Both of 
these give roughly equal contributions to the enhancement (Fig. [6]) . One could have anticipated the KK enhancement 
on simple physical grounds by arguing as follows: In Ref. 37] it was shown that the gravitational force of attraction 
between two bodies in the RS one brane model is 



gravitational force ^ < o , q' t„, q ' (63) 



That is, the Newtonian force becomes 5-dimensional on small scales. This implies that the gravitational force is 
stronger than usual on scales r < £. It then follows that the self-gravity of perturbative modes that enter the horizon 
at high energies Hi > I will be greater than those which enter at low energies H£ < 1. Therefore, we should expect 
that the amplitude of small-scale modes to be magnifled over the amplitude of large-scale modes during horizon 
crossing in braneworld cosmology, which is exactly what we have seen in our simulations.^ 



However, we should note that this simple argument neglects the influence of the KK anisotropic stress, which we know is non-negligible 
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The amplitude enhancement of perturbations is important on comoving scales < 10 AU, which are far too small 
to be relevant to present-day/cosmic microwave background measurements of the matter power spectrum. However, 
it may have an important bearing on the formation of compact objects such as primordial black holes and boson 
stars at very high energies, i.e. the greater gravitational force of attraction in the early universe will create more 
of these objects than in GR (different aspects of primordial black holes in RS cosmology in the context of various 
effective theories have been considered in Refs. [U [i^, S, S, E, S, S^l ) • We know that the abundance of primordial 
black holes can be constrained by big bang nucleosythesis and observations of high-energy cosmic rays, so it would be 
interesting to see if the kind of enhancement of density perturbations predicted in this paper can be used to derive 
new limits on Randall- Sundrum cosmology. We leave this issue for future work. 
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